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Abstract 



A dilute gas of particles with short range interactions is considered in 
a shearing stationary state. A Gaussian thermostat keeps the total kinetic 
energy constant. For infinitely many particles it is shown that the thermostat 
becomes a friction force with constant friction coefficient. For finite number 
of particles A'^, the fluctuations around this constant are of order l/^/N, and 
distributed approximately Gaussian with deviations for large fluctuations. 
These deviations prohibit a derivation of fluctuation-dissipation relations far 
from equilibrium, based on the Fluctuation Theorem. 
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The interest in the relation between non-equihbrium statistical mechanics and micro- 
scopic equations of motion, which already occupied Boltzmann, has revived in recent years, 
on the one hand due to the development of chaos theory, but even more due to results 
from non-equilibrium molecular dynamics [|l],0. The main focus in the field is on stationary 
states. A stationary state, if it is not the equilibrium state, is the result of an external 
driving force. But this force performs work on the system, so it heats up (viscous heating, 
Ohmic heating). In simulations this is often remedied by the introduction of a mechanical 
thermostat: one adds a friction force, —avi, in the equation of motion for the velocity Vi 
of each particle i, to keep the energy constant. For the thermostat variable a there are 
several choices. One could take it constant, but then one only gets a constant energy on 
average. It is also possible to have a time dependent, such that the total kinetic energy 
is constant (iso-kinetic Gaussian thermostat) or the total energy is constant (iso-energetic 
Gaussian thermostat) 0. Neither of these thermostats are very realistic, as the dissipation 
of the heating would more likely occur at the boundary, where the system is in contact 
with a heat bath, say. Other boundary formulations where the driving force and the ther- 
mostat are combined have also been studied 0,^. One hopes the choice of the thermostat 
doesn't matter in the thermodynamic limit. The equivalence of a constant a thermostat, 
the iso-kinetic thermostat and iso-energetic thermostat was proposed by Gallavotti 

The extra term in the equations of motion destroys the Liouvillian character of the flow: 
a given volume in phase space will not retain that volume. As the available phase space 
is usually finite, this means that on average over time the volume either stays constant 
(conservative case) or contracts (dissipative case). In a dissipative system a stationary state 
can exist only on a course grained scale, the dissipation continues forever but on ever finer 
scales. This dissipation happens at a rate called the phase space contraction rate which 
is proportional to the average of the thermostat variable a. This rate can be identified 
with the irreversible entropy production ^j. If we make this identification with a 
physical quantity, the precise implementation (iso-kinetic, iso-energetic, constant a, . . .) of 
the thermostat should not influence the average of a in the thermodynamic limit. Cohen 
suggested that a mechanical and a physical thermostat may give the same results as long 
as the rate of heating is much less than the rate at which heat can be transported to the 
wall and absorbed there. This suggests that when the rate of heating becomes to large, the 
thermostat does make a difference. At that point one also expects the assumption of local 
equilibrium underlying non-equilibrium thermodynamics to break down, and the entropy 
production may no longer have the form that was used to identify it with a. 

We want to know which thermostat to use for analytic treatment of dilute gases in non- 
equilibrium stationary states. As we are interested in the limit of many particles, having to 
use an a dependent on all these particles would certainly make work more difficult. In other 



analytic work on non-equilibrium stationary states, one simply takes a constant a ||10| , p 
A sketch of a proof of the equivalence of a Gaussian iso-kinetic thermostat, a Gaussian 
iso-energetic thermostat and a Nose-Hoover thermostat, was already given by Evans and 
Sarman |T^. In this paper, we will demonstrate the equivalence of an iso-kinetic Gaussian 
thermostat and a constant thermostat in the thermodynamic limit using kinetic theory on 
the Boltzmann level (i.e. at low densities) for a sheared system of short range interacting 
particles. 



2 



I. SHEARED GAS WITH SLLOD 



We consider a dilute gas of particles in d dimensions, under shear: the gas is contained 
between two plates a distance 2L apart (Fig. p. The two plates are moving in opposite 
directions with a velocity 7L. For not too large 7 one expects that a linear velocity profile 
will develop, so that the fluid velocity at y is 'yyx. 

We are interested in bulk properties, so we let L go to infinity, while the shear rate 7 
and the density p are fixed. To show the equivalence of the constant a thermostat and 
the Gaussian thermostat, it would in principle suffice just to take the horizontal dimensions 
infinite, but this way we can also move the boundary conditions to infinity. In the real 
physical system, as L get larger, the laminar flow becomes unstable and the system eventually 
develops turbulence. However, the thermostats we will consider assume the stability of the 
laminar flow [0], and suppress turbulence. 

There is a well known and often used set of equations for molecular simulations that 
describe shear, the SLLOD equations 0: 

Qi=Pi + lyd (1) 

Pi = Fi- jpiyX - api, (2) 

in which qi,Pi are the phase space coordinates of particle i, a is the thermostat variable 
and Fi represents the forces between the particles. The mass of the particles is taken to be 
one. The shear rate is constant, but a is not. It is constructed such that the kinetic energy 
K = J2i l|Pj||^/2 in the system is exactly constant, which gives 

1 ^ ^ 

This is the iso- kinetic Gaussian thermostat. Note that a depends on the positions and 
momenta of all the particles. 

The interpretation of Eq. (|l|) is that pi is the peculiar velocity of particle i with respect 
to the local fluid element that has velocity jyx. In the laboratory frame, a reasonable set 
of equations to write down is 

= Vi, Vi = Fi-a {vi - -fyix) . 

This particular form of the thermostat term is chosen because a linear velocity proflle is 
expected, and we want the temperature to be constant, so the kinetic energy in the frame 
that moves with that local velocity is to be dissipated. A Boltzmann equation for the one- 
particle distribution function will give an appropriate description at low densities. This 
equation has to be supplemented by the boundary condition that the average velocities at 
the boundaries y = ±L are ±7L. To get rid of the L dependence, one can transform the 
velocities to peculiar velocities: Pi = Vi — 'yyiX. The average (peculiar) velocity now has 
to be zero at the boundaries, so when L ^ 00 they have to be zero at inflnity: this is the 
same boundary condition as for the standard Boltzmann equation. The transformation to 
peculiar velocities yields the SLLOD equations: 

qi = Vi= Pi + 'jyiX, 
Pi = Vi- -fijiX = Fi- api - ^PiyX. 
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II. KINETIC APPROACH 



A. Effective motion of the thermostat 



The Gaussian thermostat involves an a which depends on the position and velocity 
of every particle, and so it varies in time. We'll now derive equations of motion for the 
thermostat for which we do not need to know all these positions and velocities, by introducing 
one extra thermostat variable /?. The derivation is in two steps: first we consider free flight, 
and then we take into account the effect of collisions. 

— * 

During free flight F is zero so the thermostat is given by 

Using the equations of motions, one flnds that the time derivative is 

« = -2«' + ^E^^l = -2«' + /5 (4) 

i=l 

where we've deflned the last part as the second thermostat variable /3: 

We combine these to the thermostat vector 6 = {a, (3). The time derivative of j3 is express- 
able again in terms of a and (3: 



There is a conserved quantity 



$ = -2af3. (6) 



H = 



After change of variables to X = 1/(2/5) and P = a/ (5, this conserved quantity takes on a 
Hamiltonian form H{X, P) = ^P"^ — X. The equations of motion are X — P and P — 1. 
The general solution, transformed back to 6, has the form 

^^'^^ -2H+\t-tori'~i'')^ 

with to a constant. 

So far, we only considered free flight. The duration of a free flight is very small in a 
system of N particles: it is of the order of {Nu)~^, where u is the coUision frequency of a 
single particle. In a collision, changes by an amount: 

^n^Zl( P'^^P'^y + ^2a;P2j, - PlxPly " P2xP2y 
2K \ -7 [P% + P2y - ply - ply 

^ i{Pi,P2,n) 
N 
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where primes denote the value of the variables after collisions and n is the collision parameter, 
i.e. h = {ff2i ~P2i)/||P2i ~P2i||- As K is of order A^, AO is of order A^"-*^. itself is of order 
one. The definition makes $, of order one. During a typical free flight time of one particle, 
every particle in the system has collided once on average. In each collision two particles are 
involved, so N/2 collisions have taken place in the system during this time. Thus, during 
one free flight time, the thermostat got N/2 changes of order A^~^ and this adds up to an 
effect of order one, because, as we will see, the average of A0 is non-zero. 

We will see later that the effect of the thermostat depending on all the particles - an un- 
physical idea in some sense - results just in a fluctuating thermostat with well defined mean 
and distribution. We are interested in this distribution. It will depend on the distribution 
of A0, which depends on the velocity distribution in the system, which in turn is affected 
by the thermostat distribution. But some general properties can already be derived without 
this subtle interplay. 

We want to write down a Boltzmann equation for the probability distribution function 
F{0; t) of the thermostat: 



dt 30 ^ ' dt 



(7) 



Here is given by Eq. (^) and Eq. (^). The collision integral counts the number of states 
that are lost and gained in collisions: 



dF 
'dt 



d0*dpidp2P{pi,P2,0*;t) J dn 



rate of {pi,p2,h) { 6 \ - 0* - \ - 6 {0 - 0* 

in which P is the joint distribution of pi, p2 and 0, and 

rate of {pi,p2,n) = —pB{n,p2i) 

where p2i = P2— Pi- B is the rate of collisions with n given that the colliding particles have 
relative velocity p2i- It can be expressed in the differential cross section s{n,p2i), which 
measures the number of deflected particles per unit solid angle around ^^21 when a beam of 
particles of unit density is incident on one other particle. In d dimensions: 

B{n,p2i) = 11^21 II s(n,p2i)2'^~-^ \n ■ P2if~^ 

where P21 = P21/IIP21 1| • The last factor is the Jacobian that arises because we integrate over n 
while s is defined per unit solid angle of p'2i- Strictly, we ought to take P2— Pi + ix{y2 — Ui) 
instead of P21, but in the Boltzmann approach the two particles are taken at the same 
position when they collide. In an Enskog-type approach this would matter. 

To proceed, we need an expression for P in terms of F and the one particle distribution 
function /(p), which we will take normalized to one. For /, we can also write down a 
Boltzmann equation, so we will have a system of two coupled Boltzmann equations for / and 
F . To derive the standard Boltzmann equation for /, one uses the Stosszahlansatz that 
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states that the two-particle distribution function f2{'Pi-,f>2) is proportional to the product of 
the one particle distribution functions f (pi) f {P2) , i-e. that the two particles are uncorrelated 
when they collide. We want a generalization of the Stosszahlansatz for P, but setting 
P(pi,P2, 0) = f{Pi)f{P2)F{d) can't quite be right for the following reason. 

The Stosszahlansatz can be generated by taking the phase space density p{{pi}) to be 
YliLifipi), i-e. all the particles arc uncorrelated (arguably, this is too strong a condition, 
but it will serve to make our point). Let 

f) - ^ ( -^PyP- \ 



such that 



1 ^ 

&-l;^T.Up^) 

i=l 



then the quantity 



n(pi,P2, 0) = j dps... dpNp{{Pi})S 



e 



N-2 



= fiPi)fiP2M0) 

factorizes, because the delta function doesn't involve pi and p2. One easily derives that 

N f / ^ Ne-es{pi)-es{p2)\ 



P{pi,P2,e) = {^^} n U,P2, 



N-2 



and this doesn't factorizc. 

P does however factorize to zeroth order when we expand in N~^. To see this, we first 
expand the expression for P in terms of the factorized H: 
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P(P1,P2,0)=/(P1)/(P2){$O 

+ 1 [2d - ds{p,) - d,{p2)] • 1^ } + ^(^ 



F{d) is given by / dpidp2P{pi,P2,6) so integrating the above formula gives a relation be- 
tween F and 

m = m + ^ [^m + 2 - (0.)] • + o{n-% 

where (dg) — f dpf{p)6s{p). This relation can be inverted up to order A^~^: 
m = F{d) - 1 |4F(0) + 2[d- (0,)] • + 0{N-'). 
When we put this back into the the formula for P, we get P expressed in terms of / and F: 
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P{PuP2,e) = f{p,)f{p2)[F{e) + 



■2\ 



This will serve as our Stosszahlansatz. We insert it in the collision integral and perform the 
6* integration: 



dF 
'dt 



{f ^1 dF 

J dpidp2dnf{pi)f{p2)pB{n,p2i)- \ ■ 

+ J^{J dpidp2dnf{pi)f{p2)pB{n,p2i)^ x 



y<mPi)+0sm-2{6s)]] (8) 
where we expanded in A^~^ once more. We define the coUisional averages: 

(&) ^^i^^20?^/(pi)/(p2)p5(n,_p2i)C(pi,P2,^). 
The Boltzmann equation for the thermostat can now be written as 

This is of the form of Eq. (|^) in which points in the 6 phase space follow the effective free 
flight dynamics 



6 



-2o? + /3 + a 



This amounts to just adding the average effect of collisions to the change of Q. 

Some typical trajectories in Q phase space are plotted in Fig. ^, both with and without 
the effect of collisions. When (a, 6) 7^ 0, it is no longer possible to derive the equations 
of motion from a Hamiltonian: the motion is dissipative, and there exists a flxed point 
^0 = (q^Oj/^o); deflned by 

- + /5o + a = ; 2ao/3o = h. (9) 

Physically, one expects the system to heat up without a thermostat, so the thermostat should 
act as a friction: is positive. To consider the stability of the flxed point, we linearize the 
equation of motion around it. Writing = + <^^, "we get 



The eigenvalues of the matrix are 3a;o± yc^Q — /Sq is positive, so the flxed point is stable 
only if tto > 0, consistent with the physical expectation. Furthermore, from the deflnitions 
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in Eq. and Eq. (||), one can see that < Po, so the eigenvalues are complex and the fixed 
point is a stable spiral. If we take an initial distribution within the domain of attraction, all 
points will end up in this fixed point. Thus, in the stationary state, the distribution of the 
thermostat is a delta function at 

F{0) = 5(0 - 0o) 



B. Boltzmann equation for the one particle distribution function 

As we mentioned, the analysis is not complete without a second Boltzmann equation, for 
the one particle distribution function f{p]t). The difficult part a priori is what to do with 
the free fiight term —api. The analysis of the thermostat variables however showed that the 
a is, in the stationary case, a constant, so we just replace a by ao'. 

df d 

'dt ^ ^^["^^^'^""o^^^ = J{fJ). 
where we considered again only the uniform case. For the stationary case we are interested 



in, also the time derivative disappears. The collision integral is [|13 

JUJ)\p=p^ = j dp2dnpB{n,p2i) x 



In some analytic work p!0|JTl| , a constant a is also used and its value is determined by 

demanding that the second moment = J WpW^ f {p) dp is a conserved quantity. This 
yields 

a = (11) 

Note that is constant and equal to 2K/N, so this choice of a is very analogous to 

Eq. @. In fact the solution of Eq. (j^) determining the thermostat fixed point, is this a, 
together with 

B = ^2 (Py) . (12) 

^ ^ {\\m ^ ' 



To show this we consider a and h: 

-7 



dpidp2dhpB{n,p2i)f{pi)f{p2) x 



2K/N 

X I {p'lxP'ly + P'2xP2y " PlxPly - P2xP2y] 

dpidp2dn pB{h, p2i)pixPiy x 



2 

-7 

^{f{Pl)f{P2)-f{Pl)f{p2)} 

-7 



dpJU^ f)PxPy, 



^dpJUJ)pl. 
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We insert for J(/, /) the left hand side of the Boltzmann equation and find after partial 
integration 

O {PxPy) 2 yl) o n 

imn imn 
imn 

Combined with Eq. we see that indeed (ao, Po) = («, P) given by Eq. ( ]TT| ) and Eq. 
is the self-consistent solution. 



III. THERMOSTAT FLUCTUATIONS 

When N is finite, there are fluctuations around the thermostat value ao- One can see that 
the first correction, on the right hand side of Eq. (P) has a diffusive form, with a diffusion 
coefficient of order A^~^. Combined with the drift towards a fixed point from the left hand 
side, this will make the thermostat distribution F peaked around this point, with a width 
of order N~^/'^ . For large the width is so small that the linearized equation, Eq. (|10|) is a 
good approximation for most of the distribution. For a linear fixed point, the distribution 
would be Gaussian. So the distribution will become Gaussian for large iV, except 

in the tails, where the linearized equation does not hold. 

The Gaussian nature of the distribution also carries over to the finite time average of a. 



a 



1 

— I a(t)dt. 

T Jo 



For large times r, the thermostat will spend a long time in the neighborhood of the fixed 
point (aoj Po), so will also be Gaussian distributed for large A^, but, again, with deviations 
for large fiuctuations. This all is in accordance with what one would expect from a central 
limit theorem. 

Bonetto et al. find a Gaussian distribution for a, but mention that the Gaussian 



was not what they expected, in fact, that it couldn't be Gaussian, because it would give a 
kind of generalized fiuctuation-dissipation relation (see also [|15[) when combined with the 



fiuctuation theorem of Gallavotti and Cohen 16 



However, the deviations from the Gaussian nature at large fiuctuations prohibit this 
derivation of far-from-equilibrium fiuctuation-dissipation relations. The fiuctuation theorem 
states that the probability 7r(p) of finding that a is pao, divided by the probability that it 
is — pao; satisfies for large r, 

exp [rdNaQp] . 



The result was found for Anosov systems, but has been found in numerical simulations [17 



and in Langevin equations [T^ too, and seems to have a broader validity |jT9|. Combined 



with a Gaussian form of vr, one gets a relation between the variance and the mean of the 
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distribution, i.e. a kind of fluctuation-dissipation relation. The variance can be linked to a 



correlation function and (generalized) Green-Kubo formulas emerge [14|. But for these to 



hold, the Gaussian character has to be guaranteed for negative a also, and the central limit 
theorem, nor an extension of the analysis given here, could justify that. Recently this was 
acknowledged by Searles and Evans and the Green-Kubo relations far from equilibrium 
were refuted in their molecular dynamics simulations. This does not totally explain the 
results of Bonetto et al |T^, who get a Gaussian ■Krip) fo^' just 10 particles, the Gaussian 



clearly covering negative p as well. There are of course uninvestigated prefactors in the 
widths of the Gaussian and in range of validity. 

Only when we are near equilibrium the fluctuation theorem and the central limit theorem 
can be combined, provided we flrst take the limit of the external fleld going to zero before 
we take the thermodynamic limit, even though we can take a large number of particles 
throughout. In such a scheme, the distribution tt{p) is centered almost around zero, such 
that it always apphes to some negative p. 



IV. CONCLUSIONS 

We have treated the Gaussian thermostat in a sheared system of short range interacting 
particles at low density using kinetic equations. We found that in the thermodynamic limit, 
the thermostat force becomes a friction force with a constant friction coefficient. The value 
of this constant was shown to be consistent with the requirement that the second moment 
of the one particle distribution function is conserved. This conclusion did not depend on 
a smallness of the shear rate, so it applies also far from equilibrium. The constant friction 
force has been used in other work ||TO| , p]T| and the results from those should apply equally to 
Gaussian thermostatted systems, in the thermodynamic limit. 

We briefly discussed flnite corrections. These give rise to fluctuations of the friction 
coefficient around its mean, of the order of l/\fN. The fluctuations are close to Gaussian 
for large N , except for very large fluctuations away from the mean. Far-from-equilibrium 
Green-Kubo relations rely in one derivation at least [|I4| on the fluctuation theorem ||16|| , 



which concerns large fluctuations. The breakdown of such Green-Kubo relations recently 
found by Searles and Evans pO| has its origin in the deviations from the Gaussian nature 



for large fluctuations. Near equilibrium one doesn't need the large fluctuations. 



ACKNOWLEDGEMENTS 

The author thanks Prof. H. van Beijeren, Prof. J.R. Dorfman and Prof. E.G.D. Cohen 
for valuable discussions. He acknowledges the hospitality and support of the Institute of 
Science and Technology at the University of Maryland in College Park and of the Rockefeller 
University in New York. The work reported here was supported by FOM, SMC and by 
the NWO Priority Program Non-Linear Systems, which are flnancially supported by the 
"Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)". 



10 



REFERENCES 



M. Mareschal (ed.), Special issue on the Proceedings of the Euroconference on The Mi- 
croscopic Approach to Complexity in Non-Equilibrium Molecular Simulations CECAM 
at ENS-Lyon, France, 15-19 July 1996, Physica A 240 (1997). 

D. J. Evans and G. P. Morriss, Statistical Mechanics of Nonequilihrium Liquids, (Aca- 
demic Press, London, 1990). 

N.I. Chernov and J.L. Lebowitz, Phys. Rev. Lett. 75, 2831-2834 (1995). 
Ch. Dellago and H. A. Posch, J. Stat. Phys. 88, 825-842 (1997). 
G. Gallavotti, Chaos 8, Issue 2, 384-392 (1998). 

D. Ruelle, J. Stat. Phys. 85, 1-23 (1996). 

E. G .DCohen, Physica A 240, 43-53 (1997). 

S. R. de Groot and P. Mazur, N on- Equilibrium Thermodynamics, (Dover Publications, 
Inc. New York, 1984). 

E. G. D. Cohen, Physica A 213, 293-314 (1995). 
M. Lee and J. W. Dufty, Phys. Rev. E 56, 1733-1745 (1997). 
J. F. Lutsko, Phys. Rev. Lett. 78, 243-246 (1997). 
D. J. Evans and S. Sarman, Phys. Rev. E 48, 65-70 (1993). 

C. Cercignani, Theory and Application of the Boltzmann Equation (Scottish Academic 
Press, Edinburgh/London, 1975). 

F. Bonetto, G. Gallavotti and P. L. Garrido, Physica D 105, 226-252 (1997). 

D. J. Evans and D. J. Searles, Phys. Rev. E 52, 5839-5848 (1995). 

G. Gallavotti and E. G. D. Cohen, J. Stat. Phys. 80, 931-970 (1995). 
D. J. Evans, E. G. D. Cohen and G. P. Morriss, Phys. Rev. Lett. 71, 2401-2404 (1993). 
J. Kurchan, J. Phys. A, Math. Gen. 31, 3719-3729 (1998). 

G. E. Crooks, The Gallavotti- Cohen Fluctuation Theorem and the Nonequilibrium 
Work Relation for Free Energy Differences, |http:/ /xxx.lanl.gov eprint archive, cond- 
I mat/990155^ (1999). 

[20] D. J. Searles and D. J. Evans, The Fluctuation Theorem and Green-Kubo Relations, 
Ihttp:// xxx.lanl.gov] eprint archive, |cond-mat /990202l| (1999). 



11 



FIGURES 



2L 







— > / 

X A 





FIG. 1. Velocity profile in a gas under shear 
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FIG. 2. Some typical trajectories of the thermostat. The first figure is without collisions, 
the second one is with collisional averages a and b set to the arbitrary values of —0.18 and 0.04, 
respectively. 
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